home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 275_02 / prob41.c < prev    next >
Encoding:
C/C++ Source or Header  |  1980-01-01  |  24.3 KB  |  873 lines
  1.  
  2. /* prob41.c                        */
  3. /* program for lcau41 option 't'            */
  4. /* calculate probabilities related to evolution        */
  5. /* Harold V. McIntosh, 10 August 1987            */
  6. /* 4 September 1987 - tetrahedral stereopair        */
  7. /* 30 December 1987 - contoured stereopair layers    */
  8.  
  9. /* references:                            */
  10. /*                                */
  11. /* W. John Wilbur, David J. Lipman and Shihab A. Shamma        */
  12. /* On the prediction of local patterns in cellular automata    */
  13. /* Physica 19D 397-410 (1986)                    */ 
  14. /*                                */
  15. /* Howard A. Gutowitz, Jonathan D. Victor and Bruce W. Knight    */
  16. /* Local structure theory for cellular automata            */
  17. /* Physica 28D 18-48 (1987)                    */
  18.  
  19. /*    Copyright (C) 1987    */
  20. /*    Copyright (C) 1988    */
  21. /*    Harold V. McIntosh    */
  22. /*    Gerardo Cisneros S.    */
  23.  
  24. # define BROW       14   /* row for bar charts        */
  25. # define EROW       23   /* row for evolution synopsis    */
  26.  
  27. /* edit the probability screen    */
  28.  
  29. edtri() {char c; int i;
  30.   videomode(COLGRAF);
  31.   videopalette(YELREGR);
  32.  
  33.   videocursor(0,1,33);
  34.   videoputc('2',3);
  35.   videocursor(0,13,23);
  36.   videoputc('0',3);
  37.   videocursor(0,13,39);
  38.   videoputc('1',3);
  39.  
  40.   tbase();
  41.   edges(50,3);
  42.  
  43.   while (0<1) {
  44.   videocursor(0,0,0);
  45.   hscrrul();
  46.   videocursor(0,0,36);
  47.   videoputc('?',2);
  48.   c=kbdin();
  49.   if (c == '\015') break;
  50.   videocursor(0,0,38);
  51.   videoputc(c,2);
  52.   videocursor(0,0,36);
  53.   videoputc(' ',0);
  54.   switch (c) {
  55.     case '+': videopalette(WHCYMAG); break;
  56.     case '-': videopalette(YELREGR); break;
  57.     case 'B': tbase(); break;
  58.     case 'E': tedge(); break;
  59.     case 'a': nblclr(); asfreq(3); break;
  60.     case 'e': pevolve(); break;
  61.     case 'g': for (i=0, asctobin(); i<50; i++, onegen(AL)) lifreq(2); break;
  62.     case 'G': for (i=0, asctobin(); i<200; i++, onegen(AL)) lifreq(1); break;
  63.     case 'h': edges(80,3); break;
  64.     case 'k': videohegrid(10,0); break;
  65.     case 'm': nblclr(); moncar(); break;
  66.     case 't': hpdiff(12,8); break;
  67.     case 'p': elayer(40,0); break;
  68.     case 'q': elayer(40,1); break;
  69.     case 'r': elayer(40,2); break;
  70.     case 's': elayer(40,3); break;
  71.     case 'u': elayer(72,8); break;
  72.     case 'v': olayer(72,8); break;
  73.     case 'U': elayer(42,8); break;
  74.     case 'V': olayer(42,8); break;
  75.     case '1': nblclr(); oneblfreq(8*BROW,300,56); break;
  76.     case '2': nblclr(); twoblfreq(8*BROW,300,56); break;
  77.     case '3': nblclr(); thrblfreq(8*BROW,300,56); break;
  78.     case '4': nblclr(); foublfreq(8*BROW,300,56); break;
  79.     case '5': nblclr(); fivblfreq(8*BROW,300,56); break;
  80.     case 'z': nblclr(); break;
  81.     case '/': videomode(COLGRAF); videopalette(YELREGR); break;
  82.     case '?': nblclr(); trmenu(); break;
  83.     default: break;
  84.     }; /* end switch */
  85.   };   /* end while  */
  86.   videopalette(WHCYMAG);
  87.   videomode(T80X25);
  88. }      /* end edtri  */
  89.  
  90. /* show menu */
  91. trmenu() {
  92.   videocursor(0,BROW,0);
  93.   printf("a     - a priori estimates\n");
  94.   printf("m,g,G - sample evolution\n");
  95.   printf("pqrs  - prob of states\n");
  96.   printf("t,u,v - sq, odd-even\n");
  97.   printf("12345 - n-block bar charts\n");
  98.   printf("+-    - change color quad\n");
  99.   printf("e     - 16 lines evolution\n");
  100.   printf("z/?   - clear panel, screen; show menu\n");
  101.   printf("<cr>  - exit\n");
  102. }
  103.  
  104. /* show sixteen lines of evolution on screen */
  105. pevolve() {int i, j;
  106.   videoscroll(EROW,0,EROW+1,40,0,0);
  107.   asctobin();
  108.   for (j=8*EROW; j<8*EROW+16; j++) {
  109.     for (i=0; i<AL; i++) videodot(j,i,arr1[i]);
  110.     onegen(AL);
  111.     };
  112. }
  113.  
  114. /* Clear a space for the n-block frequencies */
  115. nblclr() {videoscroll(BROW,0,BROW+8,40,0,0);}
  116.  
  117. /* locations of the weighted masses */
  118. /* for tetrahedron sitting on base  */
  119. tbase() {
  120.   wmul[0]=160.0; wmvl[0]=100.0;
  121.   wmul[1]=310.0; wmvl[1]=100.0;
  122.   wmul[2]=240.0; wmvl[2]=195.0;
  123.   wmul[3]=250.0; wmvl[3]=157.0;
  124.  
  125.   wmur[0]= 10.0; wmvr[0]=100.0;
  126.   wmur[1]=160.0; wmvr[1]=100.0;
  127.   wmur[2]= 90.0; wmvr[2]=195.0;
  128.   wmur[3]= 80.0; wmvr[3]=157.0;
  129. }
  130.  
  131. /* locations of the weighted masses */
  132. /* for tetrahedron standing on edge */
  133. tedge() {
  134.   wmul[0]=160.0; wmvl[0]=157.0;
  135.   wmul[1]=310.0; wmvl[1]=157.0;
  136.   wmul[2]=240.0; wmvl[2]=195.0;
  137.   wmul[3]=250.0; wmvl[3]=100.0;
  138.  
  139.   wmur[0]= 10.0; wmvr[0]=157.0;
  140.   wmur[1]=160.0; wmvr[1]=157.0;
  141.   wmur[2]= 90.0; wmvr[2]=195.0;
  142.   wmur[3]= 80.0; wmvr[3]=100.0;
  143. }
  144.  
  145. /* plot a single point on a triangular grid */
  146. videotrdot(u,v,x,y,z,l) double x, y, z; int u, v, l; {double s;
  147. s=199.0/(x+y+z); x=x*s; y=y*s; z=z*s;
  148. videodot(u-(int)(0.433*z),v+(int)(0.250*(y+y+z)),l);
  149. }
  150.  
  151. /* plot a single point on a tetrahedral grid */
  152. videohedot(a,b,p,l) double p[KK]; int a, b, l; {double s, u, v, w;
  153. s=199.0/(p[0]+p[1]+p[2]+p[3]);
  154. u=(p[1]+p[3])*s; v=p[2]*s; w=(p[1]-p[3])*s;
  155. videodot(a-(int)(0.433*p[1]*s),b-50+(int)(0.250*(u+u+v)),l);
  156. videodot(a-(int)(0.433*p[1]*s),b+50+(int)(0.250*(w+w+v)),l);
  157. }
  158.  
  159. /* plot a stereo point on a triangular grid */
  160. /* (u,v) screen location lower left corner  */
  161. /* (x,y,z) triangular coordinates        */
  162. /* t       height                */
  163. /* l       color                */
  164. videostdot(u,v,x,y,z,t,l) double x, y, z, t; int u, v, l; {double s;
  165. s=199.0/(x+y+z); x=x*s; y=y*s; z=z*s;
  166. videodot(u-(int)(0.433*z),v+(int)(0.250*(y+y+z+t)),l);
  167. }
  168.  
  169. /* set up a triangular gridwork */
  170. /* (u,v) = corner of triangle    */
  171. /* n = number of rulings    */
  172. /* l = color of lines        */
  173. videotrgrid(u,v,n,l) int u, v, n, l; {
  174. int i, j;
  175. double a, b, c;
  176.  
  177. if(n==0) return;
  178. for (i=0; i<=n;  i++) {
  179. for (j=0; i+j<=n; j++) {
  180.   a=((double)(i)/(double)(n));
  181.   b=((double)(j)/(double)(n));
  182.   c=1.0-a-b;
  183.   videotrdot(u,v,a,b,c,l);
  184.   };};
  185. }
  186.  
  187. /* display state frequencies in arr1 as a dot in Bezier diagram */
  188. lifreq(l) int l; {int i, stat[KK]; double staf[KK], s;
  189. s=1.0/((double)(AL));
  190. for (i=0; i<KK; i++) stat[i]=0;
  191. for (i=0; i<AL; i++) (stat[arr1[i]])++;
  192. for (i=0; i<KK; i++) staf[i]=s*((double)(stat[i]));
  193. videomassdot(staf,l);
  194. }
  195.  
  196. /* plot a single point as a weighted mass center */
  197. videomassdot(x,l) double *x; int l; {double ul, vl, ur, vr; int i;
  198. if (l>=KK) return;
  199. for (i=0, ul=0.0; i<KK; i++) ul+=x[i]*wmul[i];
  200. for (i=0, vl=0.0; i<KK; i++) vl+=x[i]*wmvl[i];
  201. for (i=0, ur=0.0; i<KK; i++) ur+=x[i]*wmur[i];
  202. for (i=0, vr=0.0; i<KK; i++) vr+=x[i]*wmvr[i];
  203. videodot(199-(int)vl,(int)ul,l);
  204. videodot(199-(int)vr,(int)ur,l);
  205. }
  206.  
  207. /* plot a point in 4 triangles - face projections of a tetrahedron */
  208. video4dots(x,l) int l; double x[KK]; {double u[KK-1];
  209. u[0]=x[0]; u[1]=x[1]; u[2]=x[2]; videotrdot(99,200,u,l);
  210. u[0]=x[1]; u[1]=x[2]; u[2]=x[3]; videotrdot(199,0,u,l);
  211. u[0]=x[2]; u[1]=x[3]; u[2]=x[0]; videotrdot(149,100,u,l);
  212. u[0]=x[3]; u[1]=x[0]; u[2]=x[1]; videotrdot(199,200,u,l);
  213. }
  214.  
  215. /* set up a triangular gridwork of edge n in color l */
  216. videohegrid(n,l) int n, l; {
  217. int i0, i1, i2;
  218. double x[KK], nn;
  219.  
  220. if (n==0) return;
  221. nn=1.0/((double)(n));
  222. for (i0=0, x[0]=0.0; i0<=n; i0++, x[0]+=nn) {
  223. for (i1=0, x[1]=0.0; i0+i1<=n; i1++, x[1]+=nn) {
  224. for (i2=0, x[2]=0.0; i0+i1+i2<=n; i2++, x[2]+=nn) {
  225.   x[3]=1.0-x[0]-x[1]-x[2];
  226.   videomassdot(x,l);
  227.   };};};
  228. }
  229.  
  230. /* show the skeleton framework of the tetrahedron */
  231. /* n = points per edge                  */
  232. /* l = color of edge                  */
  233. edges(n,l) int n, l; {int i; double x[KK], nn;
  234. if (n==0) return;
  235. nn=1.0/((double)(n));
  236. for (i=0, uvec(x,0); i<=n; i++, x[0]-=nn, x[1]+=nn) videomassdot(x,l);
  237. for (i=0, uvec(x,0); i<=n; i++, x[0]-=nn, x[2]+=nn) videomassdot(x,l);
  238. for (i=0, uvec(x,0); i<=n; i++, x[0]-=nn, x[3]+=nn) videomassdot(x,l);
  239. for (i=0, uvec(x,1); i<=n; i++, x[1]-=nn, x[2]+=nn) videomassdot(x,l);
  240. for (i=0, uvec(x,1); i<=n; i++, x[1]-=nn, x[3]+=nn) videomassdot(x,l);
  241. for (i=0, uvec(x,2); i<=n; i++, x[2]-=nn, x[3]+=nn) videomassdot(x,l);
  242. }
  243.  
  244. /* unit vector */
  245. uvec(x,i) double *x; int i; {int j;
  246. for (j=0; j<KK; j++) x[j]=0.0;
  247. x[i]=1.0;
  248. }
  249.  
  250. /* homogeneous vector */
  251. hvec(x) double *x; {int i; double s;
  252. s=1.0/((double)(KK));
  253. for (i=0; i<KK; i++) x[i]=s;
  254. }
  255.  
  256. /* layered contour map in stereoscopic tetrahedron */
  257. /* n = length of edge  */
  258. /* l = function option */
  259. elayer(n,l) int n, l; {int m, nn;
  260.   m=n/4; nn=4*m;
  261.   pdiff(3*m,nn,l);
  262.   pdiff(2*m,nn,l);
  263.   pdiff(1*m,nn,l);
  264.   pdiff(0*m,nn,l);
  265.   edges(nn,3);
  266. }
  267.  
  268. /* layered contour map in stereoscopic tetrahedron */
  269. /* n = length of edge  */
  270. /* l = function option */
  271. olayer(n,l) int n, l; {int m, nn;
  272.   m=n/8; nn=8*m;
  273.   pdiff(7*m,nn,l);
  274.   pdiff(5*m,nn,l);
  275.   pdiff(3*m,nn,l);
  276.   pdiff(1*m,nn,l);
  277.   edges(nn,3);
  278. }
  279.  
  280. /* "Monte Carlo" determination of probabilities */
  281. moncar() {
  282. int    i, b[KK];
  283. double bf[KK];
  284.   asctobin();
  285.   onegen(AL);
  286.   for (i=0; i<KK; i++) b[i]=0;
  287.   for (i=0; i<AL; i++) b[arr1[i]]++;
  288.   for (i=0; i<KK; i++) bf[i]=((double)(b[i]))/((double)(AL));
  289.   videocursor(0,BROW+8,0);
  290.   printf("(m-c) "); 
  291.   for (i=0; i<KK; i++) printf("%2d:%5.3f",i,bf[i]);
  292. }
  293.  
  294. /* Generate coefficients of the Bernstein Polynomial */
  295. berncoef() {int i, i0, i1, i2;
  296. for ( i=0;  i<KK;  i++) {
  297. for (i0=0; i0<KK; i0++) {
  298. for (i1=0; i1<KK; i1++) {
  299. for (i2=0; i2<KK; i2++) {
  300.   bp[i][i0][i1][i2]=0.0;
  301.   };};};};
  302. for (i0=0; i0<KK; i0++) {
  303. for (i1=0; i1<KK; i1++) {
  304. for (i2=0; i2<KK; i2++) {
  305.   bp[ascrule[i0][i1][i2]-'0'][i0][i1][i2]+=1.0;
  306.   };};};
  307. }
  308.  
  309. /* evaluate the nth generation Bernstein polynomial at point p */
  310. double bern(i,p) int i; double p[KK]; {int i0, i1, i2; double s;
  311. s=0.0;
  312. for (i0=0; i0<KK; i0++) {
  313. for (i1=0; i1<KK; i1++) {
  314. for (i2=0; i2<KK; i2++) {
  315.   s+=p[i0]*p[i1]*p[i2]*bp[i][i0][i1][i2];
  316.   };};};
  317. return s;
  318. }
  319.  
  320. /* graph the probability differential for state l over a triangle */
  321. pdiff(i0,n,l) int i0, n, l; {
  322. int    i1, i2;
  323. double p[KK], en, bern(), sqsu();
  324. if (n==0) return;
  325. if (i0>n) return;
  326. en=1.0/((double)(n));
  327. p[0]=((double)(i0))*en;
  328. berncoef();
  329. for (i1=0, p[1]=0.0; i0+i1<=n; i1++, p[1]+=en) {
  330. for (i2=0, p[2]=0.0; i0+i1+i2<=n; i2++, p[2]+=en) {
  331.   p[3]=1.0-p[0]-p[1]-p[2];
  332.   switch(l) {
  333.     case 0: videomassdot(p,psign(bern(0,p)-p[0])); break;
  334.     case 1: videomassdot(p,psign(bern(1,p)-p[1])); break;
  335.     case 2: videomassdot(p,psign(bern(2,p)-p[2])); break;
  336.     case 3: videomassdot(p,psign(bern(3,p)-p[3])); break;
  337.     case 4: videomassdot(p,lico(bern(0,p),i1+i2)); break;
  338.     case 5: videomassdot(p,lico(bern(1,p),i2)); break;
  339.     case 6: videomassdot(p,lico(bern(2,p),i1)); break;
  340.     case 7: videomassdot(p,lico(bern(3,p),i1)); break;
  341.     case 8: videomassdot(p,sqco(sqsu(p),i1+i2)); break;
  342.     default: break;
  343.     };
  344.   if (kbdst()) {kbdin(); return;};
  345.   };};
  346. }
  347.  
  348. /* graph the probability differential for state l over a tetrahedron    */
  349. /* n = number of points per edge                    */
  350. /* l = color of dot                            */
  351.  
  352. hpdiff(n,l) int n, l; {
  353. int    i, j, k;
  354. double bern(), x[KK], nn;
  355.  
  356. if (n==0) return;
  357. nn=1.0/((double)(n));
  358. berncoef();
  359. for (i=0, x[0]=0.0; i<=n; i++, x[0]+=nn) {
  360. for (j=0, x[1]=0.0; i+j<=n; j++, x[1]+=nn) {
  361. for (k=0, x[2]=0.0; i+j+k<=n; k++, x[2]+=nn) {
  362.   x[3]=1.0-x[0]-x[1]-x[2];
  363.   switch(l) {
  364.     case 0: videomassdot(x,psign(bern(0,x)-x[0])); break;
  365.     case 1: videomassdot(x,psign(bern(1,x)-x[1])); break;
  366.     case 2: videomassdot(x,psign(bern(2,x)-x[2])); break;
  367.     case 3: videomassdot(x,psign(bern(3,x)-x[3])); break;
  368.     case 4: videomassdot(x,lico(bern(0,x),i+j)); break;
  369.     case 5: videomassdot(x,lico(bern(1,x),j)); break;
  370.     case 6: videomassdot(x,lico(bern(2,x),i)); break;
  371.     case 7: videomassdot(x,lico(bern(3,x),i)); break;
  372.     case 8: videomassdot(x,sqco(sqsu(x),i+j)); break;
  373.     default: break;
  374.     };
  375.   if (kbdst()) {kbdin(); return;};
  376.   };};};
  377. }
  378.  
  379. /* choose a color for the contour level based on sign: red, yellow, green */
  380. psign(x) double x; {if (x>0.01) return 1; if (x<-0.01) return 2; return 3;}
  381.  
  382. /* calculate the squaresum of the four polynomials at a point */
  383. double sqsu(p) double p[KK]; {int i; double s, t, bern();
  384. for (i=0, s=0.0; i<KK; i++) {t=bern(i,p)-p[i]; s+=t*t;};
  385. return s;
  386. }
  387.  
  388. /* generate squaresum compatible contour values */
  389. int sqco(t,i) int i; double t; {int l;
  390. if (t<0.0075) l=0; else {
  391. if (t<0.0750) {if (i%2) l=0; else l=3;} else {
  392. if (t<0.1250) l=3; else {
  393. if (t<0.2500) {if (i%2) l=3; else l=2;} else {
  394. if (t<0.5000) l=2; else {
  395. if (t<0.7500) {if (i%2) l=2; else l=1;} else {
  396. if (t<1.0000) l=1; else {
  397.           if (i%2) l=1; else l=0;
  398.    }}}}}}}
  399. return l;
  400. }
  401.  
  402. /* generate linear contour values */
  403. int lico (t,i) int i; double t; {int l;
  404. if (t<0.1) l=0; else {
  405. if (t<0.2) {if (i%2) l=0; else l=3;} else {
  406. if (t<0.3) l=3; else {
  407. if (t<0.4) {if (i%2) l=3; else l=2;} else {
  408. if (t<0.5) l=2; else {
  409. if (t<0.6) {if (i%2) l=2; else l=1;} else {
  410. if (t<0.7) l=1; else {
  411.           if (i%2) l=1; else l=0;
  412.    }}}}}}}
  413. return l;
  414. }
  415.  
  416. /* graph the frequencies of ascrule in color l */
  417. asfreq(l) int l; {
  418. int    i, j, i0, i1, i2;
  419. int    stat[KK], stal, stac, star;        /* statistic counts */
  420. double staf[KK], pp;
  421.  
  422. pp=1.0/((double)(KK*KK*KK));
  423. stal=0; stac=0; star=0;
  424. for (i=0; i<KK; i++) stat[i]=0;
  425. for (i0=0; i0<KK; i0++) {
  426. for (i1=0; i1<KK; i1++) {
  427. for (i2=0; i2<KK; i2++) {
  428.   j=ascrule[i0][i1][i2]-'0';
  429.   stat[j]++;
  430.   if(j==i0) star++;
  431.   if(j==i1) stac++;
  432.   if(j==i2) stal++;
  433.   };};};
  434. videocursor(0,BROW,0);
  435. for (i=0; i<KK; i++) printf("%2d    - %5.2f\n",i,((double)stat[i])*pp);
  436. printf("\n");
  437. printf("left  - %5.2f\n",((double)stal)*pp);
  438. printf("still - %5.2f\n",((double)stac)*pp);
  439. printf("right - %5.2f\n",((double)star)*pp);
  440. for (i=0; i<KK; i++) staf[i]=((double)stat[i])*pp;
  441. videomassdot(staf,l);
  442. }
  443.  
  444. /* evaluate the one-block probabilities after one generation */
  445. onebl(x,a) double x[KK], a[KK]; {int i0, i1, i2, j0;
  446.  
  447. for (j0=0; j0<KK; j0++) x[j0]=0.0;
  448.  
  449. for (i0=0; i0<KK; i0++) {
  450. for (i1=0; i1<KK; i1++) {
  451. for (i2=0; i2<KK; i2++) {
  452.   j0=ascrule[i0][i1][i2]-'0';
  453.   x[j0]+=a[i0]*a[i1]*a[i2];
  454.   };};};
  455. }
  456.  
  457. /* iterate the 1-block parameters to find self-consistent values */
  458. /* graph the iterative steps in bar-chart form             */
  459. /* ll - initial line    */
  460. /* mm - length of line    */
  461. /* nn - number of lines    */
  462.  
  463. oneblfreq(ll,mm,nn) int ll, mm, nn; {
  464. int    ii, i, l, m, n;
  465. double op[KK], np[KK];
  466. double d, e, f, s;
  467.  
  468. m=0;
  469. f=(double)mm;
  470. s=1.0/((double)(KK));
  471. n=(int)(f*s);
  472. videodot(ll,m++,3);
  473. for (i=0; i<KK; i++) {op[i]=s; for (l=0; l<n; l++) videodot(ll,m++,i);};
  474. videodot(ll,m++,3);
  475. videomassdot(op,2);
  476.  
  477. for (ii=1; ii<=nn; ii++) {
  478.   e=0.0; m=0;
  479.   onebl(np,op);
  480.   videodot(ll+ii,m++,3);
  481.   for (i=0; i<KK; i++) {
  482.     n=(int)(f*np[i]); if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i);
  483.     d=op[i]-np[i];
  484.     e+=d*d;
  485.     op[i]=np[i];
  486.     };
  487.   videodot(ll+ii,m++,3);
  488.   if (op[0]<=0.001) break;
  489.   if (op[1]<=0.001) break;
  490.   if (op[2]<=0.001) break;
  491.   if (op[3]<=0.001) break;
  492.   if (e<=0.0000001) break;
  493.   if (kbdst()) {kbdin(); break;};
  494.   videomassdot(op,1);
  495.   };
  496. videocursor(0,BROW+8,0);
  497. printf("(1-blk)"); 
  498. for (i=0; i<KK; i++) printf("%2d:%5.3f",i,op[i]);
  499. videomassdot(op,3);
  500. }
  501.  
  502. /* evaluate the two-block probabilities after one generation */
  503. twobl(x,a) double x[KK][KK], a[KK][KK]; {
  504. int    i0, i1, i2, i3;
  505. int    j0, j1;
  506. double w, b[KK];
  507.  
  508. for (j0=0; j0<KK; j0++) {for (j1=0; j1<KK; j1++) x[j0][j1]=0.0;};
  509.  
  510. for (j0=0; j0<KK; j0++) {b[j0]=0.0; for (j1=0; j1<KK; j1++) b[j0]+=a[j0][j1];};
  511.  
  512. for (i0=0; i0<KK; i0++) {
  513. for (i1=0; i1<KK; i1++) {
  514. for (i2=0; i2<KK; i2++) {
  515. for (i3=0; i3<KK; i3++) {
  516.   j0=ascrule[i0][i1][i2]-'0';
  517.   j1=ascrule[i1][i2][i3]-'0';
  518.   w=a[i0][i1]*a[i1][i2]*a[i2][i3];
  519.   if (w!=0.0) w/=b[i1]*b[i2];
  520.   x[j0][j1]+=w;
  521.   };};};};
  522. }
  523.  
  524. /* iterate the 2-block parameters to find self-consistent values */
  525. /* graph the iterative steps in bar-chart form             */
  526. /* ll - initial line    */
  527. /* mm - length of line    */
  528. /* nn - number of lines    */
  529. twoblfreq(ll,mm,nn) int ll, mm, nn; {
  530. int    ii, i, j, l, m, n;
  531. double op[KK][KK], np[KK][KK];
  532. double b[KK], d, e, f, s;
  533.  
  534. m=0;
  535. f=(double)mm;
  536. s=1.0/((double)(KK*KK));
  537. n=(int)(f*s);
  538. videodot(ll,m++,3);
  539. for (i=0; i<KK; i++) {
  540. for (j=0; j<KK; j++) {
  541.   op[i][j]=s;
  542.   for (l=0; l<n; l++) videodot(ll,m++,j);
  543.   };};
  544. videodot(ll,m++,3);
  545. hvec(b); videomassdot(b,2);
  546.  
  547. for (ii=1; ii<=nn; ii++) {
  548.   e=0.0; m=0;
  549.   twobl(np,op);
  550.   videodot(ll+ii,m++,3);
  551.   for (i=0; i<KK; i++) {
  552.   for (j=0; j<KK; j++) {
  553.     n=(int)(f*np[i][j]);
  554.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,j);
  555.     d=op[i][j]-np[i][j];
  556.     e+=d<0.0?-d:d;
  557.     op[i][j]=np[i][j];
  558.     };};
  559.   videodot(ll+ii,m++,3);
  560. for (i=0; i<KK; i++) {b[i]=0.0; for (j=0; j<KK; j++) b[i]+=op[i][j];};
  561. videomassdot(b,1);
  562.   if (e<=0.0001) break;
  563.   if (kbdst()) {kbdin(); break;};
  564.   };
  565. videocursor(0,BROW+8,0);
  566. printf("(2-blk)"); 
  567. for (i=0; i<KK; i++) printf("%2d:%5.3f",i,b[i]);
  568. videomassdot(b,3);
  569. }
  570.  
  571. /* evaluate the three-block probabilities after one generation */
  572. thrbl(x,a) double x[KK][KK][KK], a[KK][KK][KK]; {
  573. int    i0, i1, i2, i3, i4;
  574. int    j0, j1, j2;
  575. double w, b[KK][KK];
  576.  
  577. for (j0=0; j0<KK; j0++) {
  578. for (j1=0; j1<KK; j1++) {
  579.   for (j2=0; j2<KK; j2++) x[j0][j1][j2]=0.0;
  580.   b[j0][j1]=0.0;
  581.   for (j2=0; j2<KK; j2++) b[j0][j1]+=a[j0][j1][j2];  
  582.   };};
  583.  
  584. for (i0=0; i0<KK; i0++) {
  585. for (i1=0; i1<KK; i1++) {
  586. for (i2=0; i2<KK; i2++) {
  587. for (i3=0; i3<KK; i3++) {
  588. for (i4=0; i4<KK; i4++) {
  589.   j0=ascrule[i0][i1][i2]-'0';
  590.   j1=ascrule[i1][i2][i3]-'0';
  591.   j2=ascrule[i2][i3][i4]-'0';
  592.   w=a[i0][i1][i2]*a[i1][i2][i3]*a[i2][i3][i4];
  593.   if (w!=0.0) w/=b[i1][i2]*b[i2][i3];
  594.   x[j0][j1][j2]+=w;
  595.   };};};};};
  596. }
  597.  
  598. /* iterate the 3-block parameters to find self-consistent values */
  599. /* ll - initial line    */
  600. /* mm - length of line    */
  601. /* nn - number of lines    */
  602. thrblfreq(ll,mm,nn) int ll, mm, nn; {
  603. int    ii, i, j, k, l, m, n;
  604. double op[KK][KK][KK], np[KK][KK][KK];
  605. double b[KK], bb[KK][KK], d, e, f, s;
  606.  
  607. m=0;
  608. f=(double)mm;
  609. s=1.0/((double)(KK*KK*KK));
  610. n=(int)(f*s);
  611. videodot(ll,m++,3);
  612. for (i=0; i<KK; i++) {
  613. for (j=0; j<KK; j++) {
  614. for (k=0; k<KK; k++) {
  615.   op[i][j][k]=s;
  616.   for (l=0; l<n; l++) videodot(ll,m++,k);
  617.   };};};
  618. videodot(ll,m++,3);
  619. hvec(b); videomassdot(b,2);
  620.  
  621. for (ii=1; ii<=nn; ii++) {
  622.   e=0.0; m=0;
  623.   thrbl(np,op);
  624.   videodot(ll+ii,m++,3);
  625.   for (i=0; i<KK; i++) {
  626.   for (j=0; j<KK; j++) {
  627.   for (k=0; k<KK; k++) {
  628.     n=(int)(f*np[i][j][k]);
  629.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,k);
  630.     d=op[i][j][k]-np[i][j][k];
  631.     e+=d<0.0?-d:d;
  632.     op[i][j][k]=np[i][j][k];
  633.     };};};
  634.   videodot(ll+ii,m++,3);
  635.  
  636.   for (i=0; i<KK; i++) {
  637.     b[i]=0.0;
  638.     for (j=0; j<KK; j++) {
  639.     for (k=0; k<KK; k++) {
  640.       b[i]+=op[i][j][k];
  641.       };}; };
  642.   videomassdot(b,1);
  643.   if (e<=0.0001) break;
  644.   if (kbdst()) {kbdin(); break;};
  645.   };
  646.  
  647. videocursor(0,BROW+8,0);
  648. printf("(3-blk)"); 
  649. for (i=0; i<KK; i++) printf("%2d:%5.3f",i,b[i]);
  650. videomassdot(b,3);
  651.  
  652. for (i=0; i<KK; i++) {
  653. for (j=0; j<KK; j++) {
  654.   bb[i][j]=0.0;
  655.   for (k=0; k<KK; k++) {
  656.     bb[i][j]+=op[i][j][k];
  657.     };}; };
  658.  
  659. /*videocursor(0,BROW+9,0);*/
  660. /*for (i=0; i<KK; i++) for (j=0; j<KK; j++)*/ 
  661. /*  printf("%1d%1d:%5.3f ",i,j,bb[i][j]);*/
  662. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  663. }
  664.  
  665. /* evaluate the four-block probabilities after one generation */
  666. foubl(x,a) double x[KK][KK][KK][KK], a[KK][KK][KK][KK]; {
  667. int    i0, i1, i2, i3, i4, i5;
  668. int    j0, j1, j2, j3;
  669. double w, b[KK][KK][KK];
  670.  
  671. for (j0=0; j0<KK; j0++) {
  672. for (j1=0; j1<KK; j1++) {
  673. for (j2=0; j2<KK; j2++) {
  674.   for (j3=0; j3<KK; j3++) x[j0][j1][j2][j3]=0.0;
  675.   b[j0][j1][j2]=0.0;
  676.   for (j3=0; j3<KK; j3++) b[j0][j1][j2]+=a[j0][j1][j2][j3];  
  677.   };};};
  678.  
  679. for (i0=0; i0<KK; i0++) {
  680. for (i1=0; i1<KK; i1++) {
  681. for (i2=0; i2<KK; i2++) {
  682. for (i3=0; i3<KK; i3++) {
  683. for (i4=0; i4<KK; i4++) {
  684. for (i5=0; i5<KK; i5++) {
  685.   j0=ascrule[i0][i1][i2]-'0';
  686.   j1=ascrule[i1][i2][i3]-'0';
  687.   j2=ascrule[i2][i3][i4]-'0';
  688.   j3=ascrule[i3][i4][i5]-'0';
  689.   w=a[i0][i1][i2][i3]*a[i1][i2][i3][i4]*a[i2][i3][i4][i5];
  690.   if (w!=0.0) w/=b[i1][i2][i3]*b[i2][i3][i4];
  691.   x[j0][j1][j2][j3]+=w;
  692.   };};};};};};
  693. }
  694.  
  695. /* iterate the 4-block parameters to find self-consistent values */
  696. /* ll - initial line    */
  697. /* mm - length of line    */
  698. /* nn - number of lines    */
  699. foublfreq(ll,mm,nn) int ll, mm, nn; {
  700. int    ii, i0, i1, i2, i3, l, m, n;
  701. double op[KK][KK][KK][KK], np[KK][KK][KK][KK];
  702. double b[KK], bb[KK][KK], d, e, f, s;
  703.  
  704. m=0;
  705. f=(double)mm;
  706. s=1.0/((double)(KK*KK*KK*KK));
  707. n=(int)(f*s);
  708. videodot(ll,m++,3);
  709. for (i0=0; i0<KK; i0++) {
  710. for (i1=0; i1<KK; i1++) {
  711. for (i2=0; i2<KK; i2++) {
  712. for (i3=0; i3<KK; i3++) {
  713.   op[i0][i1][i2][i3]=s;
  714.   for (l=0; l<n; l++) videodot(ll,m++,i3);
  715.   };};};};
  716. videodot(ll,m++,3);
  717. hvec(b); videomassdot(b,2);
  718.  
  719. for (ii=1; ii<=nn; ii++) {
  720.   e=0.0; m=0;
  721.   foubl(np,op);
  722.   videodot(ll+ii,m++,3);
  723.   for (i0=0; i0<KK; i0++) {
  724.   for (i1=0; i1<KK; i1++) {
  725.   for (i2=0; i2<KK; i2++) {
  726.   for (i3=0; i3<KK; i3++) {
  727.     n=(int)(f*np[i0][i1][i2][i3]);
  728.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i3);
  729.     d=op[i0][i1][i2][i3]-np[i0][i1][i2][i3];
  730.     e+=d<0.0?-d:d;
  731.     op[i0][i1][i2][i3]=np[i0][i1][i2][i3];
  732.     };};};};
  733.   videodot(ll+ii,m++,3);
  734.  
  735.   for (i0=0; i0<KK; i0++) {
  736.     b[i0]=0.0;
  737.     for (i1=0; i1<KK; i1++) {
  738.     for (i2=0; i2<KK; i2++) {
  739.     for (i3=0; i3<KK; i3++) {
  740.       b[i0]+=op[i0][i1][i2][i3];
  741.       };};}; };
  742.   videomassdot(b,1);
  743.   if (e<=0.0001) break;
  744.   if (kbdst()) {kbdin(); break;};
  745.   };
  746.  
  747. videocursor(0,BROW+8,0);
  748. printf("(4-blk)"); 
  749. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f",i0,b[i0]);
  750. videomassdot(b,3);
  751.  
  752. for (i0=0; i0<KK; i0++) {
  753. for (i1=0; i1<KK; i1++) {
  754.   bb[i0][i1]=0.0;
  755.   for (i2=0; i2<KK; i2++) {
  756.   for (i3=0; i3<KK; i3++) {
  757.     bb[i0][i1]+=op[i0][i1][i2][i3];
  758.     };}; };};
  759.  
  760. /*videocursor(0,BROW+9,0);*/
  761. /*for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++)*/ 
  762. /*  printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);*/
  763. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  764. }
  765.  
  766. /* evaluate the five-block probabilities after one generation */
  767. fivbl(x,a) double x[KK][KK][KK][KK][KK], a[KK][KK][KK][KK][KK]; {
  768. int    i0, i1, i2, i3, i4, i5, i6;
  769. int    j0, j1, j2, j3, j4;
  770. double w, b[KK][KK][KK][KK];
  771.  
  772. for (j0=0; j0<KK; j0++) {
  773. for (j1=0; j1<KK; j1++) {
  774. for (j2=0; j2<KK; j2++) {
  775. for (j3=0; j3<KK; j3++) {
  776.   for (j4=0; j4<KK; j4++) x[j0][j1][j2][j3][j4]=0.0;
  777.   b[j0][j1][j2][j3]=0.0;
  778.   for (j4=0; j4<KK; j4++) b[j0][j1][j2][j3]+=a[j0][j1][j2][j3][j4];  
  779.   };};};};
  780.  
  781. for (i0=0; i0<KK; i0++) {
  782. for (i1=0; i1<KK; i1++) {
  783. for (i2=0; i2<KK; i2++) {
  784. for (i3=0; i3<KK; i3++) {
  785. for (i4=0; i4<KK; i4++) {
  786. for (i5=0; i5<KK; i5++) {
  787. for (i6=0; i6<KK; i6++) {
  788.   j0=ascrule[i0][i1][i2]-'0';
  789.   j1=ascrule[i1][i2][i3]-'0';
  790.   j2=ascrule[i2][i3][i4]-'0';
  791.   j3=ascrule[i3][i4][i5]-'0';
  792.   j4=ascrule[i4][i5][i6]-'0';
  793.   w=a[i0][i1][i2][i3][i4]*a[i1][i2][i3][i4][i5]*a[i2][i3][i4][i5][i6];
  794.   if (w!=0.0) w/=b[i1][i2][i3][i4]*b[i2][i3][i4][i5];
  795.   x[j0][j1][j2][j3][j4]+=w;
  796.   };};};};};};};
  797. }
  798.  
  799. /* iterate the 5-block parameters to find self-consistent values */
  800. /* ll - initial line    */
  801. /* mm - length of line    */
  802. /* nn - number of lines    */
  803. fivblfreq(ll,mm,nn) int ll, mm, nn; {
  804. int    ii, i0, i1, i2, i3, i4, l, m, n;
  805. double op[KK][KK][KK][KK][KK], np[KK][KK][KK][KK][KK];
  806. double b[KK], bb[KK][KK], d, e, f, s;
  807.  
  808. m=0;
  809. f=(double)mm;
  810. s=1.0/((double)(KK*KK*KK*KK*KK));
  811. n=(int)(f*s);
  812. videodot(ll,m++,3);
  813. for (i0=0; i0<KK; i0++) {
  814. for (i1=0; i1<KK; i1++) {
  815. for (i2=0; i2<KK; i2++) {
  816. for (i3=0; i3<KK; i3++) {
  817. for (i4=0; i4<KK; i4++) {
  818.   op[i0][i1][i2][i3][i4]=s;
  819.   for (l=0; l<n; l++) videodot(ll,m++,i4);
  820.   };};};};};
  821. videodot(ll,m++,3);
  822. hvec(b); videomassdot(b,2);
  823.  
  824. for (ii=1; ii<=nn; ii++) {
  825.   e=0.0; m=0;
  826.   fivbl(np,op);
  827.   videodot(ll+ii,m++,3);
  828.   for (i0=0; i0<KK; i0++) {
  829.   for (i1=0; i1<KK; i1++) {
  830.   for (i2=0; i2<KK; i2++) {
  831.   for (i3=0; i3<KK; i3++) {
  832.   for (i4=0; i4<KK; i4++) {
  833.     n=(int)(f*np[i0][i1][i2][i3][i4]);
  834.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i4);
  835.     d=op[i0][i1][i2][i3][i4]-np[i0][i1][i2][i3][i4];
  836.     e+=d<0.0?-d:d;
  837.     op[i0][i1][i2][i3][i4]=np[i0][i1][i2][i3][i4];
  838.     };};};};};
  839.   videodot(ll+ii,m++,3);
  840.   for (i0=0; i0<KK; i0++) {
  841.     b[i0]=0.0;
  842.     for (i1=0; i1<KK; i1++) {
  843.     for (i2=0; i2<KK; i2++) {
  844.     for (i3=0; i3<KK; i3++) {
  845.     for (i4=0; i4<KK; i4++) {
  846.     b[i0]+=op[i0][i1][i2][i3][i4];
  847.     };};};}; };
  848.   videomassdot(b,1);
  849.   if (e<=0.0001) break;
  850.   if (kbdst()) {kbdin(); break;};
  851.   };
  852. videocursor(0,BROW+8,0);
  853. printf("(5-blk)"); 
  854. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f",i0,b[i0]);
  855. videomassdot(b,3);
  856.  
  857. for (i0=0; i0<KK; i0++) {
  858. for (i1=0; i1<KK; i1++) {
  859.   bb[i0][i1]=0.0;
  860.   for (i2=0; i2<KK; i2++) {
  861.   for (i3=0; i3<KK; i3++) {
  862.   for (i4=0; i4<KK; i4++) {
  863.     bb[i0][i1]+=op[i0][i1][i2][i3][i4];
  864.     };}; };};};
  865.  
  866. /*videocursor(0,BROW+9,0);*/
  867. /*for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++) */
  868. /*  printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);*/
  869. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  870. }
  871.  
  872. /* end */
  873.